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Abstract 

A class of explicit multistage time-stepping schemes is used to construct 
an algorithm for solving the compressible Navier-Stokes equations. Flexibility 
in treating arbitrary geometries is obtained with a finite-volume formulation. 
Numerical efficiency is achieved by employing techniques for accelerating 
convergence to steady state. Computer processing is enhanced through 
vectorization of the algorithm. The scheme is evaluated by solving laminar 
and turbulent flows over a flat plate and an NACA 0012 airfoil. Numerical 
results are compared with theoretical solutions or other numerical solutions 
and/or experimental data. 
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Introduction 


In recent years, much progress has been made in the development of faster 
methods for solving the compressible Navier-Stokes equations. This progress 
is due primarily to improved numerical algorithms and the arrival of 
specialized computers such as vector processors. However, many more 
contributions in these areas will be necessary before efficient and reliable 
solutions of the Navier-Stokes equations can be obtained for complex 
aerodynamic configurations . 

The current availability of a variety of solvers for the Euler equations 
provides many opportunities for exploring possible Navier-Stokes solvers. One 
strong candidate for a rapid solver is the very efficient explicit multistage 
time-stepping scheme for the Euler equations developed by Jameson, Schmidt, 
and Turkel.^ This algorithm is formally second-order accurate except near 
shock waves where the controlled addition of dissipation permits shock 
capturing without oscillations. It has the highly desirable property that the 
steady-state solution is independent of the time step. The efficiency and 
robustness of this finite-volume scheme has been demonstrated by several 
investigators Also, this scheme has been successfully applied to three- 
dimensional flow problems. ^ In the present work, this algorithm is extended 
to allow the computation of viscous flows. Computational efficiency is 
achieved for high Reynolds number viscous flows in three ways: 1) use of 
local time stepping; 2) extension of local stability range by implicit 
residual smoothing; 3) vectorization of computer code. 

While the current effort was in progress, Agarwal and Deese^ presented a 
Runge-Kutta scheme for the thin-layer Navier-Stokes equations. Their 
principal purpose was to demonstrate the versatility of a Runge-Kutta scheme 



by applying it to a variety of fluid dynamics problems. Thus, the focus of 
their work is quite different from that of the present investigation, where 
the emphasis is on accuracy and convergence acceleration. 

In this paper the governing equations and basic elements of the viscous 
flow solver are presented. An initial evaluation of the numerical method is 
achieved by solving laminar and turbulent flows over a flat plate and an NACA 
0012 airfoil. Computer results are compared to theoretical solutions or other 
numerical solutions and/or experimental data. Convergence behavior of the 
scheme is also discussed. 


Governing Equations 

Let p, (u,v), p, E, and H denote the density, Cartesian velocity 
components, pressure, total internal energy, and total enthalpy, 
respectively. The unsteady, two-dimensional Navier-Stokes equations, 
neglecting body forces and heat sources, can be written in integral form as 
follows: 
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and e , e are unit vectors of the Cartesian coordinate system (x,y), and 
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n is a unit vector normal to the surface S enclosing the volume V • In this 
paper the working fluid is air, and it is assumed to be thermally and 
calorlcally perfect. That is, the equation of state is 


p = pRT 


( 2 ) 


where R = C - C , and the specific heats C , C are constant. The 

p v* p v 

quantities p and X are the first and second coefficients of viscosity, 

2 

respectively, and X is taken to be - -j p (Stokes hypothesis). A simple 
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power law is used to determine p. The coefficient of thermal conductivity 
(k) is evaluated using the constant Prandtl number assumption. 


Finite-Volume Formulation 

In two dimensions the volume has unit depth, and thus Eq. (1) can be 
written as 
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where Si is the region of interest, and 3G is the boundary curve. Let the 
second-order tensor H be defined by 


H = Fe + Ge 
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Then, using Cartesian coordinates Eq. (3) becomes 


|- // W dx dy + / (Fdy - Gdx) = 0. (5) 

dt a an 

The computational region is partitioned with quadrilaterals, and Eq. (5) is 
applied to each quadrilateral. This process is equivalent to performing a 
mass, momentum, and energy balance on each cell. By decoupling the temporal 
and spatial terms, a systems of ordinary differential equations is obtained. 
These equations can be solved with a variety of time-stepping schemes. 

Applying Eq. (3) to an arbitrary quadrilateral (i.e., ABCD in Fig. 1) and 
approximating the line integrals with the midpoint rule we obtain 


3t ( s ijV + LW » ' ° 


ij 


where L is a spatial discretization operator, and 
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tw ii - h ab + “bc + h cd + "da- 


The components of are now cell-averaged quantities, S^j is the area of 

the cell ABCD, and the indices (i,j) identify the cell. The vectors 
H ab , H^, ^CD* ^DA re P resent t ^ ie fluxes through the sides. For example, 


F BC Ay BC G BC ^ 


BC 


(7) 


where F^q and are the mean values of F and G on the side BC, and 



4y BC ■ y C • y B 
4x BC ‘ X C - V 

Consider the x-momentum equation. The flux associated with face BC is 
H BC ■ K,j + l/2 4y BC - '' ljj+ I /2 4* BC )(pu> 1>J+1/2 

+ P i,j+l/2 Ay BC + ^BC^vis (8) 

where 

^ H BC^vis “ ^ °x ^ BC Ay BC “ ^ X yx^BC Ax BC 
and letting 5 be any cell face quantity, 

? i,j+l/ 2 = l( 5 i,j + ? i,j+l) * 

The first-order derivatives of the viscous terms, (H BC ) vis , are evaluated with 
Green's theorem. For example, 

( Vu*i/ 2 = 'Vbc " IJ u x dxdy 

" g“ / udy, 
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where 


'Vi,!* 1 /,' ( Vbc = <6 II V xdy 


= “ IT* I udx, 
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and 9J2' is the curve A'B'C'D' (see Fig. 1). If the streamwise-like 
differences associated with the viscous flux quantities are neglected (thin- 
layer Navier-Stokes assumption), the viscous terms in the flux H nr can be 


given by 


^ H BC^ vis 


’4 '(* 
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* P A (v i, j+ l " V i,j^ * (10) 

where Ay CB = -Ay^, Ax CB = -Ax^, and P A = j + U 1>j+1 )• Note that 

with the thin-layer assumption, there are viscous contributions to fluxes 
Hg£ and Hq A only. In the present work the thin-layer assumption is 
applied. For additional details on the finite-volume formulation just 


described, see Ref. 8. 
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Dissipative Terms 

In order to suppress any odd-even point decoupling in the numerical 
solution and to prevent oscillations in the vicinity of shock waves and 
stagnation points, artificial dissipation terms are added to the finite-volume 
scheme. Therefore, Eq. (6) is replaced by 

(SW) + L W - VW = 0 (11) 

where V is the artificial dissipation operator, and the indices have been 
suppressed for convenience. Through numerical experiments Jameson^ 
established that an effective form for Vw is a blend of second and fourth 
differences with coefficients that depend on the local pressure gradient. 
This form is constructed in the following way: 


Pw = P W + v w 

x y 


( 12 ) 


where V ^ W and P W are the contributions associated with the two 
coordinate directions, and in conservation form 


P W = 
x 


d±+ l /2 >j d i~ V 2 , j 
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The terms on the right side of Eqs. (13a) and (13b) have a similar form; for 
example, taking At as the Courant-Friedrichs-Lewy time step, 
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The coefficients e^ 2 ^ and e^ 4 ^ are adapted to the flow and are defined as 


follows ; 
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where typical values of the constants k^ 2 ^ and k^ 4 ^ are 1/4 and 1/256, 
respectively. In the current work pH is used instead of pE in the 
dissipative terras in the energy equation. This is done so that constant H 
can be a solution of the energy equation. 

In smooth regions of the flow field, the dissipative terms are third- 
order. The dissipation becomes first-order in the neighborhood of a shock 
wave. However, this does not compromise the global second-order accuracy of 
the finite-volume scheme. Note that the fourth difference is eliminated near 
shocks since it can have a destabilizing effect on a numerical calculation. ^ 

In the case of viscous flows significant gradients in the velocity field 
exist in the region adjacent to a solid boundary. Therefore, even though the 
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artificial viscosity coefficients (k^, ) are small, the artificial 

dissipation can produce viscous-like effects of the same order as the physical 

ones. In the present work numerical experiments for turbulent flow over a 

( 2 ) 

flat plate showed significant artificial viscous effects, even with k = 0 

(4) 

and small values of k By setting the normal contribution to the 

artificial dissipation ( V W) to zero in at least the law-of-the-wall region 
of the turbulent boundary layer accurate solutions were obtained (see Results 
and Discussion section). 


Time-Stepping Scheme 

A member of a class of four-stage schemes is employed to advance the 
solution of Eq. (11) in time. This scheme takes the following form at time 
level n: 


W 
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where on the (q+l)st stage 

RW^ q) = i (LW (q) 

J 

and 

a 1 = 1/4, (*2 = 1/3, a 3 = 1/2. 

For efficiency both the physical and artificial viscous terms are evaluated at 
the first stage and frozen for the remaining stages. The scheme of Eq. (18) 
is a modified version of the classical fourth-order Runge-Kutta scheme. The 
advantage of the modified scheme is that it requires less storage of array 
quantities in computer processing, an important consideration for three- 
dimensional calculations. According to linear stability analysis on a model 
wave equation the modified algorithm has the same stability (i.e., same 
amplification factor) as the classical scheme.-* The method of Eq. (18) is 
fourth-order accurate in time only for linear equations; it is second-order 
accurate for nonlinear equations. However, since the primary objective here 
is to compute steady-state solutions, this is not viewed as a deficiency of 
the scheme. 

The Runge-Kutta time-stepping scheme is more efficient than explicit 
schemes such as unsplit MacCormack and leapfrog.* Also, this scheme has the 
desirable property that if RW n = 0, then and thus at the final 

stage of the scheme W n+1 = W n . Therefore, the steady-state solution in 
independent of the time step, and the scheme is amenable to a variety of 
techniques for accelerating steady-state convergence. 
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Acceleration Techniques 

Three methods are employed to accelerate convergence of the basic time- 
stepping scheme. These techniques are as follows: 1) local time stepping; 2) 
enthalpy damping; 3) residual smoothing. They are discussed in the subsequent 
subsections. 


Local Time Stepping 

With this technique the solution is advanced in time with a time step 
dictated by the local stability limit. In the present work the local At is 
based on the Courant-Friedrichs-Lewy (CFL) stability limit. Local time 
stepping allows faster signal propagation, and thus faster convergence. 


Enthalpy Damping 

In Ref. 1 enthalpy damping was introduced for the Euler equations. This 
method accelerates convergence through artificial damping terms. These terms 
are added to each equation in general and depend on the deviation of the local 
total enthalpy (H) from the steady-state value (Ha,). Since the total 
enthalpy is constant throughout a steady, inviscid flow field in which uniform 
free-stream conditions exist, such forcing terms are zero in the steady state. 

Under certain assumptions the total enthalpy may be taken as constant for 
viscous flows. If the dominant viscous terms are retained (boundary-layer 
type approximation), the following form of the energy equation for laminar 
flow can be derived: 



( 19 ) 
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Then, if there is no heat transfer at solid surfaces, and the Prandtl number 
(Pr) is unity, constant H is a steady-state solution of Eq. (19). In the 
case of turbulent flow a mean flow energy equation of the same form as Eq. 
(19) can be obtained if the eddy viscosity hypothesis is used (see section on 
turbulence modeling), and small work terms are neglected. Then, if the 
surfaces are adiabatic, and both the laminar and turbulent Prandtl numbers are 
unity, constant total enthalpy is a solution of the energy equation. 

Using enthalpy damping computations exhibit a more monotone convergence 
to steady state. This property is especially useful if an acceleration 
procedure such as residual smoothing is applied with the Runge-Kutta scheme. 


Residual Smoothing 

Residual smoothing was first introduced by Lerat^® for use with the Lax- 
Wendroff scheme. Jameson^ later introduced a similar technique in conjunction 
with the Runge-Kutta schemes. With this technique the stability range of the 
basic time-stepping scheme is extended. Linear stability analysis indicates 
that the Runge-Kutta scheme with residual smoothing is unconditionally stable 
if the smoothing parameter is sufficiently large. ^ However, as revealed in 
the subsequent analysis, the fastest convergence to steady state is not 
realized with a very large time step. 
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Consider the following two-step scheme: 


u (1) = u n - aAtLu 11 

( 20 ) 

(2) n . (1) 

u= u - AtLu v ' 


where L is a spatial discretization operator. The residual smoothing for 
the second stage is given by 


1 - 



n+1 


n\ (2) 

- u 1 = u 


n 


- u 


( 21 ) 


where 6^ is the standard central difference operator, and the product is 
over the number of space dimensions. Now, consider the model problem 


3 

u„ + u + eAx u = 0. 
t x xxxx 


( 22 ) 


If e = 0, a Fourier transform of Eq. (21) yields 

1 + 6 sin 2 (G - 1) 

= — iX sin 6 — aX 2 sin 2 0 


( 23 ) 


where 



0 = k Ax 
x 


an< i is a wave number. Equation (23) can be written as 
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G = 


iX sin 8 + ctX 2 sin 2 8 
1 + 0 sin 2 | 


(2A) 


Then, taking the smoothing parameter 



Eq. (24) for large 


X reduces 


to 


G( 0) » 1 


a sin 

a . 
sin 



(25) 


Thus, without an artificial viscosity the highest frequency, 9 = it, is not 
damped. 

If the artificial third-order dissipation (e * 0) is added, then Eq. 
(24) is replaced by 


G(0) = 1 + — x , Ad . + *A) — 
1 + aX 2 sin 2 j 


(26) 


where 


A = -i sin 9 + v sin 


4 0 


and v is proportional to e. As X gets large, 

v 2 

G(tt) ~ 1 + — > 1. 


(27) 


Therefore, the scheme is not stable for large X. However, if the same 
artificial viscosity is used at both stages, then Eq. (24) is replaced by 


G( 9) = 1 


iX sin 6(1 + aXA) - Xv sin^ ~ 
1 + oX 2 sin 2 -| 


(28) 


As X approaches infinity, the coefficient of the artificial viscosity goes 
to zero, and so Eq . (25) is recovered. When the artificial viscosity term is 
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frozen, the scheme remains unconditionally stable. However, when X is 
large, the highest frequency is damped only a small amount proportional to 
1/X. Therefore, the strategy is to minimize G*G rather than choose X 
large. 

The analysis just given for a two-stage time-stepping scheme suggests 
that for multistage schemes residual smoothing is needed only at alternate 
stages (i.e., after stages 2 and 4 of a four-stage scheme). This frequency of 
smoothing has been used in solving the Euler equations. However, based on 
stability considerations for parabolic problems, ^ the residual smoothing is 
applied at every stage in the present scheme for the Navier-Stokes equations. 


Boundary Conditions 

A typical domain for a viscous airfoil calculation is considered in 
describing the boundary conditions. At the body surface the boundary 
conditions are as follows: 1) no slip (velocity components u and v are 

zero); 2) adiabatic surface condition (3T/3y = 0). The reduced normal 

momentum equation 3p/3y =0 is used to compute the surface pressure. At 
the outer boundary (ABC in Fig. 2) both inflow and outflow can occur. 
Subsonic inviscid characteristic theory indicates that three quantities should 
be specified at an inflow boundary. In the case of inflow the total enthalpy 
(H), entropy (s), and tangential velocity (U) are assigned free-stream 
values. A one-dimensional flow analysis is applied normal to the boundary, 
and the normal velocity (V) is then computed by extrapolating the Riemann 
invariant 


( 29 ) 
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frora the interior of the domain. The quantity a is the speed of sound, and 
y is the ratio of specific heats. In the case of outflow the quantities 
H, s, and U are extrapolated from the interior, and the Riemann invariant is 


V 


00 


2a co 

Y-l • 


(30) 


Two methods have been investigated for treating the downstream (outflow) 
boundary AC. They are as follows: 1) static pressure is specified, (p = p^) , 

and density and velocity components (u,v) are extrapolated from the 
Interior; 2) dependent variables (W) extrapolated from the interior. In the 
airfoil cases considered, method 2 required less steps for a converged 
solution than method 1, and the solution in the vicinity of the airfoil was 
essentially the same. However, method 2 is not consistent with characteristic 
analysis. At the present time nonref lecting-type boundary conditions are 
being investigated. 

For supersonic flows free-stream conditions are specified at inflow 
boundaries, and flow quantities are extrapolated from the interior at outflow 
boundaries . 


Turbulence Modeling 

1 2 

Using mass-averaged variables the time-averaged flow equations have the 
same form as their laminar flow counterparts, except that the stress tensor is 
augmented by the Reynolds stress tensor, the heat flux vector is augmented by 
the additional heat flux terms associated with the turbulence, and additional 
mean energy dissipation terms appear (in many cases these terras can be 



-18- 


neglected). To close the time-averaged equations for turbulent flow the eddy 
hypothesis (Reynolds stress and heat flux terms are related to mean flow-field 
gradients) is used. Moreover, the eddy viscosity, which is added to the 
molecular viscosity to obtain the effective viscosity, is computed with the 
two-layer algebraic model of Baldwin and Lomax. J The Reynolds heat flux 
terras are approximated using the constant Prandtl number assumption. Thus, 


P 




(31a) 


k = 



(31b) 


where the subscripts i and t refer to laminar and turbulent. In the 
present work the Prandtl numbers are assumed to be 1.0. 


Results and Discussion 

The numerical scheme described in the preceding sections has been applied 
to the following test problems: 1) flow over a flat plate; 2) flow over an 

NACA 0012 airfoil. For each test problem, solutions were obtained for both 
laminar and turbulent flow. In the laminar cases the free-stream Mach number 
(M^) was 0.5, and the free-stream Reynolds number (Re^) was 5 x 10^. In 
the turbulent cases the range of M^ was 0.5 - 0.756, and the range of Re^ 
was 10^ - 4.01 x 10^. The reference length for the airfoil cases was the 


chord C 
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Flat Plate Boundary-Layer Flows 

For the flat plate calculations the nondimensional distance between the 
leading edge (X/L = 0) and the downstream boundary of the rectangular 
computational domain was 1. The upper boundary of the computational domain 

was at 10 6, where 6 was the boundary-layer thickness at the downstream 

boundary. In the laminar flow case the grid consisted of 60 x 40 cells. The 
cell spacing in streamwise direction on the plate was 0.05, and the minimum 
cell height in the normal direction was 10 • The grid points were 

distributed according to a geometric progression in the normal direction. A 
60 x 60 cell grid was used in the turbulent case (same spacing as laminar case 
in streamwise direction). The normal distance from the plate to the first 
cell center was 5 xlO ^ L, which corresponded to a y + value of about 2 
(y + = yu^/v, u T is the local friction velocity and v is the kinematic 
viscosity coefficient). 

Results for laminar flow over a flat plate are presented in Figs. 3a and 

3b. The present computed solution is compared with the incompressible 

boundary-layer solution of Blasius. Figure 3a shows the variation of skin 
friction (Cf) with local Reynolds number (Re x ). Note that no attempt has been 
made to resolve the leading edge region of the plate. This probably accounts 
for the differences in the skin friction distributions. Although not shown, 
the predicted velocity profiles are essentially self-similar beyond the 
initial region (roughly 0<X/L<0.225) of the flat plate. There is reasonably 
good agreement between a representative velocity profile (X/L = 0.875) and the 
Blasius result (Fig. 3b). 

In Figs. 4a - 4c results for turbulent flow over a flat plate are 


presented. The variation of C f with Re x is displayed in Fig. 4a. The 
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nuraerical prediction exhibits fairly good agreement with the empirical curve 
fit (M^ *• 0) once fully developed turbulent flow is realized. Figure 4b 
compares solutions for the velocity profile at X/L = 0.875, which were 
calculated with the following: 1) 1/7-th-power law; 2) Steger Navier-Stokes 
code^ (Beam and Warming implicit scheme); 3) present Runge-Kutta code. Both 
of the numerical results were obtained with the Baldwin-Lomax eddy viscosity 
model, and they show excellent agreement. The same grid was used in the 
calculations. Figure 4c compares the computed law of the wall and a curve fit 
for the experimental data of several investigators^. 


NACA 0012 Airfoil Flows 

A typical grid for an airfoil calculation is displayed in Fig. 5. This 

C-type grid was generated with an algebraic procedure developed by Le 

Balleur.^ In the direction approximately perpendicular to the airfoil, the 

coordinate lines are parabolas intersecting the airfoil and outer boundary at 

predetermined locations. The second family of coordinate lines were 

determined by distributing points along the parabolas in accordance with a 

prescribed function. Additional information concerning grid generation can be 

obtained from Ref. 16. For all "C" grids used herein the outer boundary was 

located 5 to 6 chords from the airfoil in all directions. In the laminar flow 

case, calculations were done on meshes with 128 x 32 cells and 128 x 64 cells. 

The cell spacing in the streamwise direction over the central part of the 

airfoil was AX/C = 0.05. The minimum cell height in the normal— like direction 

-4 -4 

for these grids was about 3 x 10 chords and 6 x 10 chords, 

respectively. For turbulent flow computations (120 x 50 cell grid) the 
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distance from the airfoil surface to the first cell center was approximately 
5 x 10 ^ chords. The chordwise spacing at the midsection of the airfoil 
was AX/C = 0.036. 

Figures 6a - 6c show results for the laminar flow (M = 0.5, 

3 

Re^ = 5 x 10 ) computation on the 128 x 64 cell grid. In Fig. 6a the 
pressure distribution for the viscous flow is compared with that obtained by 
solving the Euler equations. The presence of strong viscous effects is 
evident over the aft section of the airfoil. Moreover, as indicated in Fig. 
6b, the flow separates at X/C “ 0.817. The size of the thin recirculation 
region is indicated in the plot of streamlines in Fig. 6c. Since the velocity 
of the flow is very small in this thin region, the accuracy of the solution 
there is of special concern. The solution on the 128 x 64 cell grid was 

compared with that on a 128 x 32 cell grid, where the mesh spacing in the 
normal-like direction is about a factor of 2 larger. These solutions are 
essentially the same except in the reverse flow region. A bubble that is 

slightly larger in longitudinal extent (approximately 7 percent larger) is 
predicted with the finer grid. 

Results for subsonic, turbulent flow over the airfoil at zero angle of 
attack are presented in Figs. 7a - 7c. Figure 7a shows a comparison of the 
calculated pressure distribution with experimental data.^ In Fig. 7b the 
predicted skin friction distribution is compared with that obtained with a 

1 O 

boundary-layer code. For the boundary-layer calculation a very fine mesh 

spacing (more than an order of magnitude smaller than that in the present 

calculation) was used at the surface, and 100 points were placed across the 
boundary layer. Since the viscous-inviscid interaction effects are small for 
this flow, a single pass computation was done. The pressure distribution from 
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the thin-layer Navier-Stokes solution was prescribed • There is fairly good 
agreement except at the leading edge, where the present solution does not have 
adequate resolution for the very thin turbulent boundary layer. In Fig. 7c 
pressure contours for this case are displayed. 

Computed pressure distributions for subsonic (M = 0.5), viscous flow 
over the airfoil at two angles of attack (1.77 degrees and 3.51 degrees) are 
compared with those for inviscid flow and with experimental data in Figs. 8 
and 9. In each case the angle of attack (a) has been corrected for wind 
tunnel wall effects. There is good agreement between the present Navier- 
Stokes results and experiment. The influence of viscosity on these solutions, 
in particular at the suction peak, is evident. 

In Figs. 10a and 10b results for supercritical flow (M^ = 0.756) and 
zero angle of attack are presented. The pressure distributions from the 
viscous and inviscid calculations show good agreement with the data. The 
supersonic region of the flow is indicated by the Mach number contours in Fig. 
10b. 


Convergence 

For all computations in this paper the flow field was initialized with 
free-stream conditions. Moreover, the calculations were started impulsively 
by suddenly introducing the body into a uniform flow and immediately enforcing 
the appropriate boundary conditions at its surface. To measure convergence 
the root mean square of the residual of the continuity equation was used. In 
Fig. 11 convergence histories are shown for the laminar airfoil case on 
a 128 x 32 cell grid. For the history labeled A the local CFL number was 
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2, and enthalpy damping was used. For the history labeled B the CFL number 
was 6, and both enthalpy damping and residual smoothing were used. These 
results indicate that satisfactory convergence for engineering applications is 
achieved in 1300 time steps, which corresponds to 2.5 minutes on the Vector 
Processing System (VPS) 32 at Langley Research Center. Figures 12a and 12b 
show residual histories for the zero angle of attack, subcritical, turbulent 
airfoil case. In Fig. 12a the convergence history, which is highly 
oscillatory in character, is for a calculation using only local time stepping 
to accelerate convergence. As seen in Fig. 12b, a computation using enthalpy 
damping and residual smoothing (CFL number of 5) exhibits much more monotonic 
convergence behavior. In this case acceptable convergence is realized in about 
1250 time steps (4 minutes on the VPS 32). Note that the present computer code 
is not optimized for the VPS 32 system. With optimization the computer 
processing times can be reduced by a factor of 2 to 3. 

The final rate of reduction of the residual for the turbulent flow cases 
is quite slow. Based on previous work with Euler equations 5 , significant 
improvement in this rate appears to be possible with a multigrid scheme. At 
the present time this is being investigated. 


Concluding Remarks 

A finite-volume scheme for numerical integration of the Euler equations 
has been extended to allow solution of the Navier-Stokes equations. The new 
procedure, which is based on a class of four stage Runge-Kutta time-stepping 
schemes, has been made numerically efficient through the following convergence 
acceleration techniques: 1) local time stepping; 2) enthalpy damping; 3) 
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residual smoothing. Also, the high degree of vectorization possible with the 
algorithm has yielded an efficient program for vector processors. The scheme 
has been evaluated by solving the thin layer form of the Navier-Stokes 
equations for laminar and turbulent flows over a flat plate and an NACA 0012 
airfoil. Numerical results have compared well with either theoretical or 
other numerical solutions and/or experimental data. 



References 


* Jameson, A., Schmidt, W. , and Turkel, E.: "Numerical Solutions of the Euler 
Equations by Finite Folume Methods Using Runge-Kutta Time-Stepping 
Schemes," AIAA Paper 81-1259, June 1981. 

o 

Schmidt, W., Jameson, A., and Whitfield, D.: "Finite Volume Solution for 
the Euler Equation for Transonic Flow Over Airfoils and Wings Including 
Viscous Effects," AIAA Paper 81-1265, June 1981. 

O 

Salas, M. D. , Jameson, A., and Melnik, R. E.: "A Comparative Study of the 

Nonuniqueness Problem of the Potential Equation," AIAA Paper 83-1888, 
July 1983. 

^Deese , J. E. and Agarwal, R. K. , "Calculation of Axisymmetric Inlet 
Flowfield Using the Euler Equations," AIAA Paper 83-1853, July 1983. 

Jameson, A. and Baker, T. J. , "Solution of the Euler Equations for Complex 
Configurations," AIAA Paper 83-1929, July 1983. 

^ Agarwal , R. K. and Deese, J. E.: "Transonic Wing— Body Calculations Using 

Euler Equations," AIAA Paper 83-0501, January 1983. 

7 Agarwal, R. K. and Deese, J. E.: "Computation of Transonic Viscous Airfoil, 

Inlet, and Wing Flowfields," AIAA Paper 84-1551, June 1984. 



-26- 


^Peyret, R. and Taylor, T. D.: Computational Methods for Fluid Flow , 

Springer Series in Computational Physics, Springer-Verlag, 1983. 

^Turkel, E.: "Acceleration to a Steady State for the Euler Equations," NASA 

CR-172398 and ICASE Report No. 84-32, July 1984. 

l^Lerat, A.: "Une classe de schemas aux 

differences implicites pour les systemes hyperboliques de lois de 
conservation," Comptes Rendus Acad. Sciences Paris , Vol. 288A, 1979. 

^Jameson, A.: "The Evolution of Computational Methods in Aerodynamics," J. 

Appl. Mech., Vol. 50, 1983. 

^Rubesin, M. W. and Rose, W. C.: "The Turbulent Mean-Flow, Reynolds-Stress , 

and Heat-Flux Equations in Mass-Averaged Dependent Variables," NASA TM 
X-62,248, March 1973. 

^Baldwin, B. S. and Lomax, H. : "Thin Layer Approximation and Algebraic Model 

for Separated Turbulent Flows," AIAA Paper 78-257, January 1978. 

^Steger, J. L.: "Implicit Finite Difference Simulation of Flow About 

Arbitrary Geometries with Application to Airfoils," AIAA Paper 77—665, 
June 1977. 

^Hinze, J. 0.: Turbulence , 2nd edition, McGraw-Hill Book Company, Inc., New 

York, 1975. 



-27- 


iD Le Balleur, J. C.: "Strong Matching Method for Computing Transonic Viscous 

Flows Including Wakes and Separations - Lifting Airfoils," La Recherche 
Aerospatiale (English edition), No. 1981-3. 

17 Thibert, J. J., Granjacques, M. , and Ohman, L. H. : NACA 0012 Airfoil, AGARD 

Advisory Report No. 138, Experimental Data Base for Computer Program 
Assessment, May 1979. 

1 Q 

Miner, E. W., Anderson, E. C., and Lewis, C. H. : "A Computer Program for 
Two-Dimensional and Axisymmetric Nonreacting Perfect Gas and Equilibrium 
Chemically Reacting Laminar, Transitional, and/or Turbulent Boundary 
Layer Flows," VPI-E-71-8, May 1971 (available as NASA CR-132601). 



Figure 1. Finite-volume mesh. 


Figure 2. Computational domain for airfoil 
calculations . 




Figure 3. Laminar flat plate flow: 

(a) comparison of skin-friction distributions; 

(b) comparison of velocity profiles (M = 0.5, Re = 5 x IQ 3 ). 
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Law Of The Wall, M = 0 
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Figure 4. Turbulent flat plate flow: 

(a) comparison of skin friction distributions; 

(b) comparison of velocity profiles; 

(c) comparison of near wall velocity profiles 
(M - 0.5, Re = 10 6 ). 
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Figure 5. Partial view of typical grid for airfoil calculations. 








(c) 


Figure 7. Turbulent flow over an NACA 0012 airfoil: 

(a) pressure distributions; (b) skin friction distributions; 
(c) pressure contours (M = 0.5, Re = 2.89 x 10^). 
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Figure 8. Pressure distributions for flow over 

an NACA 0012 airfoil at angle of attack 

(M « * °* 5 > Re o> 88 2.91xi0 6 , a = 1.77°). 

Lift coefficient: Experiment - 0.195; 

Viscous - 0.198; Inviscid - 0.252. 
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Figure 9. Pressure distributions for flow over 

an NACA 0012 airfoil at angle of attack 

< M oo~ °- 5 » 2.93xl0 6 , a = 3.51°). 

Lift coefficient: Experiment - 0.386; 

Viscous - 0.388; Inviscid - 0.429. 
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Figure 10. Supercritical flow over an NACA 0012 airfoil: 

(a) Pressure distributions; (b) Mach number contours 
(M - 0.756, Re - 4.01 x 10 6 ). 
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Iterations 

Figure 11. Convergence histories for an NACA 0012 airfoil computation 
(M^ = 0.5, Re^ = 5000, 128 x 32 cell grid). 

A - enthalpy damping; B - enthalpy damping and residual smoothing. 




(a) (b) 

Figure 12. Convergence histories for NACA 0012 airfoil computation 
(M^ - 0.5, Re — - 2.89 x 10 6 , 120 x 50 cell grid): 

(a) local time-stepping only; (b) enthalpy damping and residual smoothing. 






1. Report No. NASA CR- 172527 2. Government Accession No. 

ICASE Report No. 84-62 

4. Title and Subtitle 

A MULTISTAGE TIME-STEPPING SCHEME FOR THE NAVIER- STOKES 


EQUATIONS 

6. Performing Organization Code 

7. Author(s) 

R. C. Swanson and Eli Turkel 

8. Performing Organization Report No. 

84-62 


10. Work Unit No. 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 


and Engineering 

Mail Stop 132C, NASA Langley Research Center 
Hampton, VA 23665 

11. Contract or Grant No. 

NAS1-17070 

13. Type of Report and Period Covered 

12. Sponsoring Agency Name and Address 

Contractor Report 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

14. Sponsoring Agency Code 

505-31-83-01 

15. Supplementary Notes 

Langley Technical Monitor: J. C. South, Jr. 

Final Report 

l — — 


16. Abstract 


3. Recipient’s Catalog No. 
5. Report Date 

February 1985 


A class of explicit multistage time-stepping schemes is used to construct an 
algorithm for solving the compressible Navier-Stokes equations. Flexibility in 
treating arbitrary geometries is obtained with a finite-volume formulation. Numerical 
efficiency is achieved by employing techniques for accelerating convergence to steady 
state. Computer processing is enhanced through vectorization of the algorithm. The 
scheme is evaluated by solving laminar and turbulent flows over a flat plate and an 
NACA 0012 airfoil. Numerical results are compared with theoretical solutions or 
other numerical solutions and/or experimental data. 


17. Key Words (Suggested by Author(s)) 

Navier-Stokes , 
finite volume, 
Runge-Kutta. 


18. Distribution Statement 

34 - Fluid Mechanics and Heat Transfer 
Unclassified - Unlimited 

19. Security Oassif. (of this report) 
Unclassified 

20. Security Classif. (of this page) 
Unclassified 

21. No. of Pages 
36 

22. Price 
A0 3 


For sale by the National Technical Information Service, Springfield, Virginia 22161 


NASA-Langley, 198 

















1 1 1 1 1 ill ill j mi ; lllllil II | j 3 1 1 1 


3 1176 005 


9 3645 



